# Replication Code for:
# Aspiration versus Apprehension: Economic Opportunities and Electoral Preferences
# Silja Häusermann, Thomas Kurer and Delia Zollinger
# British Journal of Political Science
# 2023


# Preparation
dev.off()
cat("\014")  
# globals
options(scipen=999)

# packages


# load glm.predict first
if (!require("glm.predict")) install.packages("glm.predict") 
library(glm.predict)

# then other packages
if (!require("pacman")) install.packages("pacman") 

pacman::p_load(tidyverse, magrittr, broom, hrbrthemes, plm, estimatr, sandwich, lmtest, AER, lfe, huxtable, margins, readstata13, texreg, reshape2, readxl, interplot, ggpubr, gridExtra, xtable, haven, corrplot, nnet, ggeffects, emmeans)

# ggplot theme
theme_set(theme_bw() + theme(text = element_text(size=14)))


# define paths

datpath <- "...set path..."
outpath <- "...set path..."
cwpath <-  "...set path..."


# IMPORT ----

# read data

op <- read.csv(paste0(datpath, "aspiration_apprehension_data.csv"))

# read crosswalk for task grousp

isco_task_cw <- read.csv(paste0(cwpath, "isco08_3d-task3.csv"))


# DATA PREPARATION ----

# vote choice last election

op$vc[op$partyfam2==6] <- "1 rad right"
op$vc[op$partyfam2==5 | op$partyfam2==4] <- "2 ms right"
op$vc[op$partyfam2==3 | op$partyfam2==2] <- "3 ms left"
op$vc[op$partyfam2==1] <- "4 rad left"

op$vcG[op$partyfam2==6] <- "1 rad right"
op$vcG[op$partyfam2==5 | op$partyfam2==4] <- "2 ms right"
op$vcG[op$partyfam2==3] <- "3 soc dem"
op$vcG[op$partyfam2==2] <- "4 green"
op$vcG[op$partyfam2==1] <- "5 rad left"

# participation in last election

op$part <- NA
op$part[op$X5_3_Participation==1] <- "yes"
op$part[op$X5_3_Participation==2] <- "no"

# incumbency

op$incumbent <- 0

op$incumbent[op$country=="Denmark" & (op$X5_4_Party_choice_last_election=="Liberal Alliance" | 
                                        op$X5_4_Party_choice_last_election=="Venstre" |
                                        op$X5_4_Party_choice_last_election=="Conservative People's Party" | 
                                        op$X5_4_Party_choice_last_election=="Danish People's Party")] <- 1

op$incumbent[op$country=="Germany" & (op$X5_4_Party_choice_last_election=="SPD" | 
                                        op$X5_4_Party_choice_last_election=="CDU/CSU")] <- 1


op$incumbent[op$country=="Ireland" & op$X5_4_Party_choice_last_election=="Fine Gael"] <- 1

op$incumbent[op$country=="Italy" & (op$X5_4_Party_choice_last_election=="LN" | 
                                        op$X5_4_Party_choice_last_election=="M5S")] <- 1

op$incumbent[op$country=="Netherlands" & (op$X5_4_Party_choice_last_election=="D66" | 
                                        op$X5_4_Party_choice_last_election=="VVD" |
                                        op$X5_4_Party_choice_last_election=="CDA" | 
                                        op$X5_4_Party_choice_last_election=="CU")] <- 1

op$incumbent[op$country=="Spain" & op$X5_4_Party_choice_last_election=="PSOE"] <- 1

op$incumbent[op$country=="Sweden" & (op$X5_4_Party_choice_last_election=="SAP" | 
                                        op$X5_4_Party_choice_last_election=="MP")] <- 1

op$incumbent[op$country=="United Kingdom" & op$X5_4_Party_choice_last_election=="Conservatives"] <- 1

op$incumbent[op$X5_4_Party_choice_last_election=="NA"] <- NA

# multinomial Y: incumbent, mainstream opposition, radical opposition
# hierarchy: mainstream vs. radical first, incumbent vs. (any) opposition second

op$y <- NA
op$y[op$vc=="2 ms right" | op$vc =="3 ms left"] <- "2 opp"
op$y[op$vc=="1 rad right" | op$vc=="4 rad left"] <- "3 rad"
op$y[op$incumbent==1 & op$y!="3 rad"] <- "1 inc" 

op$y[op$country=="Italy" & op$X5_4_Party_choice_last_election=="M5S"] <- "3 rad"

op$y_inc <- ifelse(op$y=="1 inc", 1, 0)
op$y_inc[is.na(op$y)] <- NA

op$y_opp <- ifelse(op$y=="2 opp", 1, 0)
op$y_opp[is.na(op$y)] <- NA

op$y_rad <- ifelse(op$y=="3 rad", 1, 0)
op$y_rad[is.na(op$y)] <- NA

# Y version with green party separated (for robustness)

op$yg <- NA
op$yg[op$vc=="2 ms right" | op$vc =="3 ms left"] <- "2 opp"
op$yg[op$vcG=="4 green"] <- "3 green"
op$yg[op$vc=="1 rad right" | op$vc=="4 rad left"] <- "4 rad"
op$yg[op$incumbent==1 & op$y!="3 rad"] <- "1 inc" 

# Y as factor, define reference

op$y2 <- as.factor(op$incumbent)
op$y2 <- relevel(op$y2, ref="1")

# dichotomous high/low income for quadrants 

table(op$income)
op$inchl <- NA
op$inchl[op$income <= 5] <- "low"
op$inchl[op$income > 5] <- "high"
table(op$inchl)

# high/medium/low income alternative (robustness)

table(op$income)
op$inchml <- NA
op$inchml[op$income <= 3] <- "low"
op$inchml[op$income > 3 & op$income < 8] <- "medium"
op$inchml[op$income >= 8] <- "high"
table(op$inchml)

# opportunity perceptions respondents for quadrants

op$dpessr <- ifelse(op$org>5, "positive", "negative")

# opportunity perceptions kids for quadrants (robustness)

op$dpessk <- ifelse(op$okg>5, "positive", "negative")

# opportunity perceptions labor market respondents for quadrants (robustness)

op$dpessrlm <- ifelse(op$orlm>5, "positive", "negative")

# education groups
op$educ_grp2[op$educ_grp==1] <- "low"
op$educ_grp2[op$educ_grp==2] <- "medium"
op$educ_grp2[op$educ_grp==3] <- "high"

op$educ_grp2 <-  ordered(op$educ_grp2, levels = c("low", "medium", "high"))

# class labels
op$class8[op$class8==1] <- "employers"
op$class8[op$class8==2] <- "small bus owners"
op$class8[op$class8==3] <- "technical profs"
op$class8[op$class8==4] <- "prod workers"
op$class8[op$class8==5] <- "managers"
op$class8[op$class8==6] <- "clerks"
op$class8[op$class8==7] <- "sociocult profs"
op$class8[op$class8==8] <- "service workers"

# age groups

op$age_grp[op$age_grp==1] <- "18-25"
op$age_grp[op$age_grp==2] <- "26-35"
op$age_grp[op$age_grp==3] <- "36-45"
op$age_grp[op$age_grp==4] <- "46-55"
op$age_grp[op$age_grp==5] <- "56-65"
op$age_grp[op$age_grp==6] <- "66+"

# task groups via isco 3 digit crosswalk

# isco finetuning for merging
op$isco3dm <- op$short_isco3d
op$isco3dm[op$isco3dm==200] <- 211
op$isco3dm[op$isco3dm==442] <- 441
op$isco3dm[op$isco3dm==230] <- 231
op$isco3dm[op$isco3dm==520] <- 521
op$isco3dm[op$isco3dm==210] <- 211
op$isco3dm[op$isco3dm==822] <- 821
op$isco3dm[op$isco3dm==530] <- 531
op$isco3dm[op$isco3dm==320] <- 321
op$isco3dm[op$isco3dm<0] <- NA

# crosswalk isco 88 3d -> task3

op <- merge(op, isco_task_cw, by.x="isco3dm", by.y="isco08_3d", all.x=TRUE)

# task groups

op$nrc <- ifelse(op$task==3, 1, 0)

# task groups label

op$taskgroup <- NA
op$taskgroup[op$task==1] <- "non-routine cognitive"
op$taskgroup[op$task==2] <- "routine"
op$taskgroup[op$task==3] <- "non-routine manual"

op$taskgroup <- factor(op$taskgroup, levels=c("non-routine manual", "routine", "non-routine cognitive"))

# quadrant respondent social

op$quadrant <- NA
op$quadrant[op$inchl=="high" & op$dpessr=="positive"] <- "comfortable"
op$quadrant[op$inchl=="high" & op$dpessr=="negative"] <- "apprehensive"
op$quadrant[op$inchl=="low" & op$dpessr=="positive"] <- "aspirational"
op$quadrant[op$inchl=="low" & op$dpessr=="negative"] <- "burdened"

op$quadrant <- factor(op$quadrant, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

# quadrant respondent labor market

op$quadrant_lm <- NA
op$quadrant_lm[op$inchl=="high" & op$dpessrlm=="positive"] <- "comfortable"
op$quadrant_lm[op$inchl=="high" & op$dpessrlm=="negative"] <- "apprehensive"
op$quadrant_lm[op$inchl=="low" & op$dpessrlm=="positive"] <- "aspirational"
op$quadrant_lm[op$inchl=="low" & op$dpessrlm=="negative"] <- "burdened"

op$quadrant_lm <- factor(op$quadrant_lm, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

# quadrant kids social

op$quadrant_kids <- NA
op$quadrant_kids[op$inchl=="high" & op$dpessk=="positive"] <- "comfortable"
op$quadrant_kids[op$inchl=="high" & op$dpessk=="negative"] <- "apprehensive"
op$quadrant_kids[op$inchl=="low" & op$dpessk=="positive"] <- "aspirational"
op$quadrant_kids[op$inchl=="low" & op$dpessk=="negative"] <- "burdened"

op$quadrant_kids <- factor(op$quadrant_kids, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

# quadrant respondent social version with 3 inc classes

op$quadrant3 <- NA
op$quadrant3[op$inchml=="high" & op$dpessr=="positive"] <- "highincpos"
op$quadrant3[op$inchml=="high" & op$dpessr=="negative"] <- "highincneg"
op$quadrant3[op$inchml=="medium" & op$dpessr=="positive"] <- "medincpos"
op$quadrant3[op$inchml=="medium" & op$dpessr=="negative"] <- "medincneg"
op$quadrant3[op$inchml=="low" & op$dpessr=="positive"] <- "lowincpos"
op$quadrant3[op$inchml=="low" & op$dpessr=="negative"] <- "lowincneg"

# additional covariates (robustness)

op$fiscal_constraint <- NA
op$fiscal_constraint[op$X1_3_Fiscal_constraint_realstc_==1 | op$X1_3_Fiscal_constraint_realstc_==2] <- 0
op$fiscal_constraint[op$X1_3_Fiscal_constraint_realstc_==3 | op$X1_3_Fiscal_constraint_realstc_==4] <- 1

op$ever_unemp_benefits <- 0
op$ever_unemp_benefits[op$X6_2_1_Experience_with_the_WS==1] <- 1

op$unemp <- 0
op$unemp[op$X8_8_7 == 1] <- 1

op$involpart <- 0
op$involpart[op$X8_11_voluntary_involuntary_wrk==1] <- 1

op$empstat_vuln <- ifelse(op$unemp==1 | op$involpart==1, 1, 0)

# export for multinomial analyses in stata (supplementary material)

op_stata <- op
names(op_stata) <- tolower(names(op_stata))
write.csv(op_stata, paste0(datpath, "op_for_stata.csv"), na="", row.names=FALSE)


################################################################################################################

# MAIN ANALYSIS ----

# Figure 1 ----

fig1 <- op %>% dplyr::select(org, orlm, okg, oklm)
fig1l <- pivot_longer(fig1, org:oklm)

fig1l$group <- NA
fig1l$group[fig1l$name=="org"] <- "Respondent, General"
fig1l$group[fig1l$name=="orlm"] <- "Respondent, Labour Market"
fig1l$group[fig1l$name=="okg"] <- "Kids, General"
fig1l$group[fig1l$name=="oklm"] <- "Kids, Labour Market"

fig1l$groupf <- factor(fig1l$group, levels=c("Respondent, General", "Respondent, Labour Market", "Kids, General", "Kids, Labour Market"))


ggplot(fig1l, aes(x = value)) +
  geom_histogram(colour="black", fill="salmon2", alpha=0.3, binwidth=1) +
  geom_vline(aes(xintercept = mean(value, na.rm = T)),
             colour = "red", linetype ="dotted", size = .8) +
  xlab("Distribution Opportunities") + ylab("Frequency") +
  facet_wrap(~ groupf)

ggsave(paste0(outpath, "BJPS_Figure1.pdf"), height=6, width=8)

# Figure 2 ----
# background only

ggplot(data=subset(op, !is.na(quadrant))) + 
  facet_wrap(~quadrant) +
  scale_color_discrete(name="") +
  labs(x="",
       y="",
       title="Socio-Democraphic Characteristics By Quadrant",
       subtitle="Qualitative Summary of Emblematic Cases Based on Multinomial Modelling") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Figure2.pdf"), height=8, width=8)


# Table 2 ----

op$y <- as.factor(op$y)
op$y <- relevel(op$y, ref="1 inc")

# opportunity respondent

fig3a_mn <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

# Table out

texreg(fig3a_mn, beside=TRUE, digits=3,
       label="tab:reg_party_org_mn", caption="Opportunity Respondent and Party Support (Multinomial, Reference: Incumbent Voting)", caption.above=TRUE,
       custom.model.names=c("Mainstream Opposition", "Radical"),
       custom.coef.map = list('org' = 'Social Opportunity',
                              'income' = 'Income',
                              'age'='Age',
                              'female' = 'Female',
                              'factor(educ)2' = 'Educ: Primary',
                              'factor(educ)3' = 'Educ: Secondary I',
                              'factor(educ)4' = 'Educ: Secondary II',
                              'factor(educ)5' = 'Educ: Post-Secondary',
                              'factor(educ)6' = 'Educ: Short Tertiary',
                              'factor(educ)7' = 'Educ: Tertiary I',
                              'factor(educ)8' = 'Educ: Tertiary II',
                              'factor(class5)2' = 'Class: Lower Service',
                              'factor(class5)3' = 'Class: Small Business',
                              'factor(class5)4' = 'Class: Skilled Worker',
                              'factor(class5)5' = 'Class: Unskilled Worker'),
       custom.note=paste("%stars. Multinomial Specification. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Table2.tex")) # Table 1 is theory table

# Figure 3 ----

# panel (a): opportunity respondent

fig3a_mn <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)


fig3a_mn <- broom::tidy(fig3a_mn,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn <- dplyr::filter(fig3a_mn, term=="org")

ref3a_mn <- c("1 inc", "org", 0, 0, 0, 0, 0, 0)
fig3a_label <- c("Mainstream Opposition", "Radical", "Incumbent (Reference)")

fig3a_plot <- rbind(fig3a_mn, ref3a_mn)
fig3a_plot <- cbind(fig3a_label, fig3a_plot)

fig3a_plot %<>% arrange(fig3a_label)
fig3a_plot[,c(4:9)]  <-  apply(fig3a_plot[,c(4:9)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3a_plot, aes(x=estimate,y=fig3a_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
   xlim((-0.25), 0.05) +
  xlab("Opportunity Respondent and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Figure3a.pdf"), width=6, height=4)


# panel (b): opportunity kids

fig3b_mn <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

fig3b_mn <- broom::tidy(fig3b_mn,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn <- dplyr::filter(fig3b_mn, term=="okg")

ref3a_mn <- c("1 inc", "okg", 0, 0, 0, 0, 0, 0)
fig3b_label <- c("Mainstream Opposition", "Radical", "Incumbent (Reference)")

fig3b_plot <- rbind(fig3b_mn, ref3a_mn)
fig3b_plot <- cbind(fig3b_label, fig3b_plot)

fig3b_plot %<>% arrange(fig3b_label)
fig3b_plot[,c(4:9)]  <-  apply(fig3b_plot[,c(4:9)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3b_plot, aes(x=estimate,y=fig3b_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
  xlim((-0.25), 0.05) +
  xlab("Opportunity Kids and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Figure3b.pdf"), width=6, height=4)

# Figure 4 ----

# panel (a): opportunity respondent

vc_inc <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)

fig_vc_inc <- broom::tidy(vc_inc,conf.int=TRUE, exponentiate=FALSE)
fig_vc_inc <- dplyr::filter(fig_vc_inc, term=="org")

fig_vc_inc_ref <- c("org", 0, 0, 0, 0, 0, 0)
fig_vc_inc_label <- c("Opposition", "Incumbent (Reference)")

fig_vc_inc <- rbind(fig_vc_inc, fig_vc_inc_ref)
fig_vc_inc <- cbind(fig_vc_inc_label, fig_vc_inc)
fig_vc_inc[,c(3:8)]  <-  apply(fig_vc_inc[,c(3:8)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig_vc_inc, aes(x=estimate,y=fig_vc_inc_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
   xlim((-0.25), 0.05) +
  xlab("Opportunity Respondent and Incumbency") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")

ggsave(paste0(outpath, "BJPS_Figure4a.pdf"), width=6, height=4)

# panel (b): opportunity kids

vc_inc2 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)

fig_vc_inc2 <- broom::tidy(vc_inc2,conf.int=TRUE)
fig_vc_inc2 <- dplyr::filter(fig_vc_inc2, term=="okg")

fig_vc_inc2_ref <- c("okg", 0, 0, 0, 0, 0, 0)
fig_vc_inc2_label <- c("Opposition" , "Incumbent (Reference)")

fig_vc_inc2 <- rbind(fig_vc_inc2, fig_vc_inc2_ref)
fig_vc_inc2 <- cbind(fig_vc_inc2_label, fig_vc_inc2)
fig_vc_inc2[,c(3:8)]  <-  apply(fig_vc_inc2[,c(3:8)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig_vc_inc2, aes(x=estimate,y=fig_vc_inc2_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
   xlim((-0.25), 0.05) +
  xlab("Opportunity Kids and Incumbency") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")

ggsave(paste0(outpath, "BJPS_Figure4b.pdf"), width=6, height=4)


# Table 3 ----

op$female <- as.factor(op$female)
op$age_grp <- as.factor(op$age_grp)

fig4_inc <- glm(y_inc~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_opp <- glm(y_opp~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_rad <- glm(y_rad~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Table out


texreg(list(fig4_inc, fig4_opp, fig4_rad), digits=3, 
       label="tab:reg_party_quadrant", caption="Opportunity Types and Party Support", caption.above=TRUE,
       custom.model.names=c("Incumbent", "Mainstream Opposition", "Radical"), 
       custom.coef.map = list('quadrantapprehensive' = 'Type: Apprehensive',
                              'quadrantaspirational' = 'Type: Aspirational',
                              'quadrantburdened' = 'Type: Burdened',
                              'age_grp26-35' = 'Age: 26-35',
                              'age_grp36-45' = 'Age: 36-45',
                              'age_grp46-55' = 'Age: 46-55',
                              'age_grp56-65' = 'Age: 56-65',
                              'age_grp66+'='Age: 66 and older',
                              'female1' = 'Female',
                              'educ_grp2.L' = 'Educ: Medium',
                              'educ_grp2.Q' = 'Educ: High',
                              'factor(class8)employers' = 'Class: Employers',
                              'factor(class8)managers' = 'Class: Managers',
                              'factor(class8)prod workers' = 'Class: Production Workers',
                              'factor(class8)service workers' = 'Class: Service Workers',
                              'factor(class8)small bus owners' = 'Class: Small Business Owners',
                              'factor(class8)sociocult profs' = 'Class: Soc. Cult. Profs',
                              'factor(class8)technical profs' = 'Class: Tech. Profs'),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Table3.tex"))


# Figure 5 ----

output_inc = predicts(fig4_inc, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_inc$party <- "incumbent"
output_opp = predicts(fig4_opp, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_opp$party <- "Mainstream Opposition"
output_rad = predicts(fig4_rad, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_rad$party <- "Radical"

output <- full_join(output_inc, output_opp)
output <- full_join(output, output_rad)

output <-  output %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_avg <- output %>% group_by(party, quadrant) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

output_avg$quadrant <-  ordered(output_avg$quadrant, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

ggplot(output_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=5) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Incumbent", "Mainstream Opposition", "Radical")), values=c(10,4,13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Figure5.pdf"), height=6, width=8)

# Figure 6 ----

op$inc5 <- NA
op$inc5[op$income==1 | op$income==2] <- 1
op$inc5[op$income==3 | op$income==4] <- 2
op$inc5[op$income==5 | op$income==6] <- 3
op$inc5[op$income==7 | op$income==8] <- 4
op$inc5[op$income==9 | op$income==10] <- 5

op$org5 <- NA
op$org5[op$org==0 | op$org==1 | op$org==2] <- 1
op$org5[op$org==3 | op$org==4] <- 2
op$org5[op$org==5 | op$org==6] <- 3
op$org5[op$org==7 | op$org==8] <- 4
op$org5[op$org==9 | op$org==10] <- 5

# Incumbent

interaction_inc <- glm(y_inc~inc5*org5+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")

interaction_inc_est <- ggplot_build(interplot(m = interaction_inc, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[1]][,c(5,3)]
interaction_inc_ci <- ggplot_build(interplot(m = interaction_inc, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[2]][,c(3,5,6)]

interaction_inc_plot <- cbind(interaction_inc_est, interaction_inc_ci)

ggplot(interaction_inc_plot, aes(x=x,y=y,color=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax)) + 
  ylab("Pr(Incumbent Vote)") + xlab("Income") +
  scale_colour_discrete(name="Opportunity Perception", labels = c("Negative", "Positive"))

interaction_inc_minmax <- interaction_inc_plot %>% filter(x==1 | x==5)

interaction_inc_minmax$xlab <- ifelse(interaction_inc_minmax$x==1, "low", "high")
interaction_inc_minmax$xlab <- ordered(interaction_inc_minmax$xlab, levels = c("low", "high"))

ggplot(interaction_inc_minmax, aes(x=xlab,y=y, shape=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax), size=.75, position = position_dodge(width = 0.1)) + 
  ylab("Pr(Incumbent Vote)") + xlab("Income") +
  ylim(0,70) +
  scale_shape_discrete(name="Opportunity \nPerception", labels = c("Negative", "Positive"))
ggsave(paste0(outpath, "BJPS_Figure6a.pdf"), width=5, height=4)

# Opposition

interaction_opp <- glm(y_opp~inc5*org5+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_opp_est <- ggplot_build(interplot(m = interaction_opp, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[1]][,c(5,3)]
interaction_opp_ci <- ggplot_build(interplot(m = interaction_opp, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[2]][,c(3,5,6)]

interaction_opp_plot <- cbind(interaction_opp_est, interaction_opp_ci)

ggplot(interaction_opp_plot, aes(x=x,y=y,color=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax)) + 
  ylab("Pr(Opposition Vote)") + xlab("Income") +
  scale_colour_discrete(name="Opportunity Perception", labels = c("Negative", "Positive"))

interaction_opp_minmax <- interaction_opp_plot %>% filter(x==1 | x==5)

interaction_opp_minmax$xlab <- ifelse(interaction_opp_minmax$x==1, "low", "high")
interaction_opp_minmax$xlab <- ordered(interaction_opp_minmax$xlab, levels = c("low", "high"))

ggplot(interaction_opp_minmax, aes(x=xlab,y=y,shape=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax), size=.75, position = position_dodge(width = 0.1)) + 
  ylab("Pr(Opposition Vote)") + xlab("Income") +
  ylim(0,70) +
  scale_shape_discrete(name="Opportunity \nPerception", labels = c("Negative", "Positive"))
ggsave(paste0(outpath, "BJPS_Figure6b.pdf"), width=5, height=4)

# Radical

interaction_rad <- glm(y_rad~inc5*org5+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_rad_est <- ggplot_build(interplot(m = interaction_rad, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[1]][,c(5,3)]
interaction_rad_ci <- ggplot_build(interplot(m = interaction_rad, var1 = "inc5", var2 = "org5", predPro = TRUE, var2_vals = c(min(op$org5, na.rm=T), max(op$org5, na.rm=T))))$data[[2]][,c(3,5,6)]

interaction_rad_plot <- cbind(interaction_rad_est, interaction_rad_ci)

ggplot(interaction_rad_plot, aes(x=x,y=y,color=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax)) + 
  ylab("Pr(Radical Vote)") + xlab("Income") +
  scale_colour_discrete(name="Opportunity Perception", labels = c("Negative", "Positive"))

interaction_rad_minmax <- interaction_rad_plot %>% filter(x==1 | x==5)

interaction_rad_minmax$xlab <- ifelse(interaction_rad_minmax$x==1, "low", "high")
interaction_rad_minmax$xlab <- ordered(interaction_rad_minmax$xlab, levels = c("low", "high"))

ggplot(interaction_rad_minmax, aes(x=xlab,y=y,shape=factor(group)))+
  geom_pointrange(aes(ymin=ymin,
                     ymax=ymax), size=.75, position = position_dodge(width = 0.1)) +
  ylab("Pr(Radical Vote)") + xlab("Income") +
  ylim(0, 70) +
  scale_shape_discrete(name="Opportunity \nPerception", labels = c("Negative", "Positive"))
ggsave(paste0(outpath, "BJPS_Figure6c.pdf"), width=5, height=4)

# Figure 7 ----

interaction_inc_cont <- glm(y_inc~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_opp_cont <- glm(y_opp~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_rad_cont <- glm(y_rad~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")

interplot(interaction_inc_cont, var1="income", var2="org") + 
  geom_hline(yintercept=0, color="red", linetype="dotted") +
  ggtitle("Incumbent Voting") + xlab("Income") + ylab("Slope of Opportunity Perception") +
  ylim((-0.15), 0.15)

ggsave(paste0(outpath, "BJPS_Figure7a.pdf"), width=5, height=4)


interplot(interaction_opp_cont, var1="income", var2="org") + 
  geom_hline(yintercept=0, color="red", linetype="dotted") +
  ggtitle("Opposition Voting") + xlab("Income") + ylab("Slope of Opportunity Perception") +
  ylim((-0.15), 0.15)

ggsave(paste0(outpath, "BJPS_Figure7b.pdf"), width=5, height=4)


interplot(interaction_rad_cont, var1="income", var2="org") + 
  geom_hline(yintercept=0, color="red", linetype="dotted") +
  ggtitle("Radical Voting") + xlab("Income") + ylab("Slope of Opportunity Perception") +
  ylim((-0.15), 0.15)

ggsave(paste0(outpath, "BJPS_Figure7c.pdf"), width=5, height=4)

# Figure 8 ----

fig4_rr <- glm(rr~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_msr <- glm(msr~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_msl <- glm(msl~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_rl <- glm(rl~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

output_rr = predicts(fig4_rr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_rr$party <- "Rad Right"
output_msr = predicts(fig4_msr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_msr$party <- "MS Right"
output_msl = predicts(fig4_msl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_msl$party <- "MS Left"
output_rl = predicts(fig4_rl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_rl$party <- "Rad Left"

output2 <- full_join(output_rr, output_msr)
output2 <- full_join(output2, output_msl)
output2 <- full_join(output2, output_rl)


output2 <-  output2 %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output2_avg <- output2 %>% group_by(party, quadrant) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output2_avg$quadrant <-  ordered(output2_avg$quadrant, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))
output2_avg$party <-  ordered(output2_avg$party, levels = c("Rad Left", "MS Left", "MS Right", "Rad Right"))

ggplot(output2_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=6) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Rad Left", "MS Left", "MS Right", "Rad Right")), values=c(10,4, 8, 13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Figure8.pdf"), height=6, width=8)

# Figure 9 ----

w <- op %>% filter(!is.na(quadrant)) %>%
                     mutate(one=1, wone=Weight_Age_sex___education) %>%
  group_by(country) %>% mutate(cntrytotal=sum(one), wcntrytotal=sum(Weight_Age_sex___education, na.rm=T)) %>%
  group_by(country, quadrant) %>%
  summarise(total=sum(one, na.rm=T), wtotal=sum(wone, na.rm=T), 
         share=total/cntrytotal, wshare=wtotal/wcntrytotal,
         perc=round(share*100, 1), wperc=round(wshare*100, 1)) %>%
  filter(row_number(total) == 1) %>%
    ungroup()


highs <- w %>% dplyr::filter(quadrant=="comfortable" | quadrant=="aspirational") %>% group_by(country) %>% mutate(order=sum(wshare)) %>% ungroup() %>% droplevels()
lows <- w %>% dplyr::filter(quadrant=="apprehensive" | quadrant=="burdened") %>% droplevels()
lows$quadrant <- factor(lows$quadrant, levels = c("burdened", "apprehensive"))

w$pos <- NA
w$pos[w$quadrant=="comfortable"] <- 0.35
w$pos[w$quadrant=="comfortable" & w$country=="Italy"] <- 0.32

w$pos[w$quadrant=="aspirational"] <- 0.1
w$pos[w$quadrant=="apprehensive"] <- -0.05
w$pos[w$quadrant=="burdened"] <- -0.25

ggplot() + geom_bar(data=highs, aes(x = reorder(country, order), y=wshare, fill=quadrant), position="stack", stat="identity") +
  geom_bar(data=lows, aes(x = country, y=-wshare, fill=quadrant), position="stack", stat="identity") +
  geom_hline(yintercept = 0, color =c("white")) +
  geom_text(data=w, aes(x=country, y=pos, label=wperc, 
                ), 
            position = position_dodge(width=1)) +
  coord_flip() + ylab("Relative Share (in %)") + xlab("") +
  scale_fill_grey(breaks=c("burdened","apprehensive","aspirational", "comfortable")) + 
  theme(legend.position = "bottom", axis.text.x = element_blank(), legend.title=element_blank())

ggsave(paste0(outpath, "BJPS_Figure9.pdf"), height=6, width=8)

# END MAIN ANALYSIS






# SUPPLEMENTARY MATERIAL ----

# Table B.1 ----

cor_op <- op %>% dplyr::select(org, orlm, okg, oklm)
names(cor_op) <- c("Resp_G", "Resp_LM", "Kids_G", "Kids_LM")

tab_b1 <- cor(cor_op, use="complete")
tab_b1 <- round(tab_b1, 2)
print(xtable(tab_b1, type = "latex"), file = paste0(outpath, "BJPS_Appendix_Table_B1.tex"))

# Table C.1 ----

fig3b_mn <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

texreg(fig3b_mn, beside=TRUE, digits=3,
       label="tab:reg_party_okg_mn", caption="Opportunity Kids and Party Support (Multinomial, Reference: Incumbent Voting)", caption.above=TRUE,
       custom.model.names=c("Mainstream Opposition", "Radical"),
       custom.coef.map = list('okg' = 'Kids Opportunity',
                              'income' = 'Income',
                              'age'='Age',
                              'female1' = 'Female',
                              'factor(educ)2' = 'Educ: Primary',
                              'factor(educ)3' = 'Educ: Secondary I',
                              'factor(educ)4' = 'Educ: Secondary II',
                              'factor(educ)5' = 'Educ: Post-Secondary',
                              'factor(educ)6' = 'Educ: Short Tertiary',
                              'factor(educ)7' = 'Educ: Tertiary I',
                              'factor(educ)8' = 'Educ: Tertiary II',
                              'factor(class5)2' = 'Class: Lower Service',
                              'factor(class5)3' = 'Class: Small Business',
                              'factor(class5)4' = 'Class: Skilled Worker',
                              'factor(class5)5' = 'Class: Unskilled Worker'),
       custom.note=paste("%stars. Multinomial Specification. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_C1.tex"))

# Table C.2 ----

vc_inc <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)
vc_inc2 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)

texreg(list(vc_inc, vc_inc2), digits=3,
       label="tab:reg_party_econvote", caption="Opportunity and Incumbent Support", caption.above=TRUE,
       custom.model.names=c("Opportunity Respondent", "Opportunity Kids"),
       custom.coef.map = list('org' = 'General Opportunity (Respondent)',
                              'okg' = 'General Opportunity (Kids)',
                              'income' = 'Income',
                              'age'='Age',
                              'female' = 'Female',
                              'factor(educ)2' = 'Educ: Primary',
                              'factor(educ)3' = 'Educ: Secondary I',
                              'factor(educ)4' = 'Educ: Secondary II',
                              'factor(educ)5' = 'Educ: Post-Secondary',
                              'factor(educ)6' = 'Educ: Short Tertiary',
                              'factor(educ)7' = 'Educ: Tertiary I',
                              'factor(educ)8' = 'Educ: Tertiary II',
                              'factor(class5)2' = 'Class: Lower Service',
                              'factor(class5)3' = 'Class: Small Business',
                              'factor(class5)4' = 'Class: Skilled Worker',
                              'factor(class5)5' = 'Class: Unskilled Worker'),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_C2.tex"))

# Table C.3 ----


fig4_mn_table <- multinom(y~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op)

texreg(list(fig4_mn_table), digits=3, beside=TRUE,
       label="tab:reg_party_quadrant_mn", caption="Opportunity Types and Party Support (Reference: Incumbent)", caption.above=TRUE,
       custom.model.names=c("Mainstream Opp", "Radical"), 
       custom.coef.map = list('quadrantapprehensive' = 'Type: Apprehensive',
                              'quadrantaspirational' = 'Type: Aspirational',
                              'quadrantburdened' = 'Type: Burdened',
                              'age_grp26-35' = 'Age: 26-35',
                              'age_grp36-45' = 'Age: 36-45',
                              'age_grp46-55' = 'Age: 46-55',
                              'age_grp56-65' = 'Age: 56-65',
                              'age_grp66+'='Age: 66 and older',
                              'female1' = 'Female',
                              'educ_grp2.L' = 'Educ: Medium',
                              'educ_grp2.Q' = 'Educ: High',
                              'factor(class8)employers' = 'Class: Employers',
                              'factor(class8)managers' = 'Class: Managers',
                              'factor(class8)prod workers' = 'Class: Production Workers',
                              'factor(class8)service workers' = 'Class: Service Workers',
                              'factor(class8)small bus owners' = 'Class: Small Business Owners',
                              'factor(class8)sociocult profs' = 'Class: Soc. Cult. Profs',
                              'factor(class8)technical profs' = 'Class: Tech. Profs'),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_C3.tex"))

# Table C.4 ----


interaction_inc_cont <- glm(y_inc~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_opp_cont <- glm(y_opp~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")
interaction_rad_cont <- glm(y_rad~income*org+age+female+educ+factor(class8)+factor(country), data=op, family="binomial")


texreg(list(interaction_inc_cont, interaction_opp_cont, interaction_rad_cont), digits=3, 
       label="tab:reg_interaction_cont", caption="Interacting Opportunity and Income (Continuous)", caption.above=TRUE,
       custom.model.names=c("Incumbent", "MS Opposition", "Radical"), 
       custom.coef.map = list('org' = 'Social Opportunity',
                              'income' = 'Income',
                              'income:org' = 'Income X Opportunity',
                              'age' = 'Age',
                              'female1' = 'Female',
                              'educ' = 'Education',
                              'factor(class8)employers' = 'Class: Employers',
                              'factor(class8)managers' = 'Class: Managers',
                              'factor(class8)prod workers' = 'Class: Production Workers',
                              'factor(class8)service workers' = 'Class: Service Workers',
                              'factor(class8)small bus owners' = 'Class: Small Business Owners',
                              'factor(class8)sociocult profs' = 'Class: Soc. Cult. Profs',
                              'factor(class8)technical profs' = 'Class: Tech. Profs'),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_C4.tex"))

# Table C.5 ----


fig4_rr <- glm(rr~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_msr <- glm(msr~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_msl <- glm(msl~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_rl <- glm(rl~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

texreg(list(fig4_rl, fig4_msl, fig4_msr, fig4_rr), digits=3, 
       label="tab:reg_partyfam_quadrant", caption="Opportunity Types and Party Support", caption.above=TRUE,
       custom.model.names=c("Rad Left", "MS Left", "MS Right", "Rad Right"), 
       custom.coef.map = list('quadrantapprehensive' = 'Type: Apprehensive',
                              'quadrantaspirational' = 'Type: Aspirational',
                              'quadrantburdened' = 'Type: Burdened',
                              'age_grp26-35' = 'Age: 26-35',
                              'age_grp36-45' = 'Age: 36-45',
                              'age_grp46-55' = 'Age: 46-55',
                              'age_grp56-65' = 'Age: 56-65',
                              'age_grp66+'='Age: 66 and older',
                              'female1' = 'Female',
                              'educ_grp2.L' = 'Educ: Medium',
                              'educ_grp2.Q' = 'Educ: High',
                              'factor(class8)employers' = 'Class: Employers',
                              'factor(class8)managers' = 'Class: Managers',
                              'factor(class8)prod workers' = 'Class: Production Workers',
                              'factor(class8)service workers' = 'Class: Service Workers',
                              'factor(class8)small bus owners' = 'Class: Small Business Owners',
                              'factor(class8)sociocult profs' = 'Class: Soc. Cult. Profs',
                              'factor(class8)technical profs' = 'Class: Tech. Profs'),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_C5.tex"))

# Figure D.1 ----

fig3a_mn <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)
fig3a_mn_rob1 <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(country), data = op)
fig3a_mn_rob2 <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(ever_unemp_benefits) + factor(country), data = op)
fig3a_mn_rob3 <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(empstat_vuln) + factor(country), data = op)
fig3a_mn_rob4 <- multinom(y ~ org + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), data = op)

fig3a_mn <- broom::tidy(fig3a_mn,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn <- dplyr::filter(fig3a_mn, term=="org")
fig3a_mn$model <- "main"

fig3a_mn_rob1 <- broom::tidy(fig3a_mn_rob1,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn_rob1 <- dplyr::filter(fig3a_mn_rob1, term=="org")
fig3a_mn_rob1$model <- "main + fiscal constraint"

fig3a_mn_rob2 <- broom::tidy(fig3a_mn_rob2,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn_rob2 <- dplyr::filter(fig3a_mn_rob2, term=="org")
fig3a_mn_rob2$model <- "main + unemp benefit"

fig3a_mn_rob3 <- broom::tidy(fig3a_mn_rob3,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn_rob3 <- dplyr::filter(fig3a_mn_rob3, term=="org")
fig3a_mn_rob3$model <- "main + vulnerability"

fig3a_mn_rob4 <- broom::tidy(fig3a_mn_rob4,conf.int=TRUE, exponentiate=FALSE)
fig3a_mn_rob4 <- dplyr::filter(fig3a_mn_rob4, term=="org")
fig3a_mn_rob4$model <- "main + all 3"

fig3a_mn_robust <- rbind(fig3a_mn, fig3a_mn_rob1, fig3a_mn_rob2, fig3a_mn_rob3, fig3a_mn_rob4)

ref <- c("1 inc", "org", 0, 0, 0, 0, 0, 0, "main")

fig3a_mn_robust <- rbind(ref, fig3a_mn_robust)

fig3a_mn_robust$label <- NA
fig3a_mn_robust$label[fig3a_mn_robust$y.level=="1 inc"] <- "Incumbent (Refernce)"
fig3a_mn_robust$label[fig3a_mn_robust$y.level=="2 opp"] <- "Mainstream Opposition"
fig3a_mn_robust$label[fig3a_mn_robust$y.level=="3 rad"] <- "Radical"


fig3a_mn_robust %<>% arrange(label)
fig3a_mn_robust[,c(3:8)]  <-  apply(fig3a_mn_robust[,c(3:8)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3a_mn_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Respondent and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted") +
  scale_color_grey()

fig3b_mn <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

fig3b_mn_rob1 <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(country), data = op)
fig3b_mn_rob2 <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(ever_unemp_benefits) + factor(country), data = op)
fig3b_mn_rob3 <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(empstat_vuln) + factor(country), data = op)
fig3b_mn_rob4 <- multinom(y ~ okg + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), data = op)

fig3b_mn <- broom::tidy(fig3b_mn,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn <- dplyr::filter(fig3b_mn, term=="okg")
fig3b_mn$model <- "main"

fig3b_mn_rob1 <- broom::tidy(fig3b_mn_rob1,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn_rob1 <- dplyr::filter(fig3b_mn_rob1, term=="okg")
fig3b_mn_rob1$model <- "main + fiscal constraint"

fig3b_mn_rob2 <- broom::tidy(fig3b_mn_rob2,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn_rob2 <- dplyr::filter(fig3b_mn_rob2, term=="okg")
fig3b_mn_rob2$model <- "main + unemp benefit"

fig3b_mn_rob3 <- broom::tidy(fig3b_mn_rob3,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn_rob3 <- dplyr::filter(fig3b_mn_rob3, term=="okg")
fig3b_mn_rob3$model <- "main + vulnerability"

fig3b_mn_rob4 <- broom::tidy(fig3b_mn_rob4,conf.int=TRUE, exponentiate=FALSE)
fig3b_mn_rob4 <- dplyr::filter(fig3b_mn_rob4, term=="okg")
fig3b_mn_rob4$model <- "main + all 3"

fig3b_mn_robust <- rbind(fig3b_mn, fig3b_mn_rob1, fig3b_mn_rob2, fig3b_mn_rob3, fig3b_mn_rob4)

ref <- c("1 inc", "okg", 0, 0, 0, 0, 0, 0, "main")

fig3b_mn_robust <- rbind(ref, fig3b_mn_robust)

fig3b_mn_robust$label <- NA
fig3b_mn_robust$label[fig3b_mn_robust$y.level=="1 inc"] <- "Incumbent (Refernce)"
fig3b_mn_robust$label[fig3b_mn_robust$y.level=="2 opp"] <- "Mainstream Opposition"
fig3b_mn_robust$label[fig3b_mn_robust$y.level=="3 rad"] <- "Radical"


fig3b_mn_robust %<>% arrange(label)
fig3b_mn_robust[,c(3:8)]  <-  apply(fig3b_mn_robust[,c(3:8)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3b_mn_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Kids and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()


vc_inc <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)

vc_inc_rob1 <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(country), family="binomial", data = op)
vc_inc_rob2 <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(ever_unemp_benefits) + factor(country), family="binomial", data = op)
vc_inc_rob3 <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(empstat_vuln) + factor(country), family="binomial", data = op)
vc_inc_rob4 <- glm(y2 ~ org + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), family="binomial", data = op)

vc_inc <- broom::tidy(vc_inc,conf.int=TRUE, exponentiate=FALSE)
vc_inc <- dplyr::filter(vc_inc, term=="org")
vc_inc$model <- "main"

vc_inc_rob1 <- broom::tidy(vc_inc_rob1,conf.int=TRUE, exponentiate=FALSE)
vc_inc_rob1 <- dplyr::filter(vc_inc_rob1, term=="org")
vc_inc_rob1$model <- "main + fiscal constraint"

vc_inc_rob2 <- broom::tidy(vc_inc_rob2,conf.int=TRUE, exponentiate=FALSE)
vc_inc_rob2 <- dplyr::filter(vc_inc_rob2, term=="org")
vc_inc_rob2$model <- "main + unemp benefit"

vc_inc_rob3 <- broom::tidy(vc_inc_rob3,conf.int=TRUE, exponentiate=FALSE)
vc_inc_rob3 <- dplyr::filter(vc_inc_rob3, term=="org")
vc_inc_rob3$model <- "main + vulnerability"

vc_inc_rob4 <- broom::tidy(vc_inc_rob4,conf.int=TRUE, exponentiate=FALSE)
vc_inc_rob4 <- dplyr::filter(vc_inc_rob4, term=="org")
vc_inc_rob4$model <- "main + all 3"

vc_inc_robust <- rbind(vc_inc, vc_inc_rob1, vc_inc_rob2, vc_inc_rob3, vc_inc_rob4)

ref <- c("org", 0, 0, 0, 0, 0, 0, "main")

vc_inc_robust <- rbind(ref, vc_inc_robust)

vc_inc_robust$label <- NA
vc_inc_robust$label[vc_inc_robust$conf.high==0] <- "Incumbent (Reference)"
vc_inc_robust$label[vc_inc_robust$conf.high!=0] <- "Opposition"


vc_inc_robust[,c(2:7)]  <-  apply(vc_inc_robust[,c(2:7)], 2, function(x) as.numeric(as.character(x)));

ggplot(vc_inc_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Respondent and Incumbency") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()


vc_inc2 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), family="binomial", data = op)

vc_inc2_rob1 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(country), family="binomial", data = op)
vc_inc2_rob2 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(ever_unemp_benefits) + factor(country), family="binomial", data = op)
vc_inc2_rob3 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(empstat_vuln) + factor(country), family="binomial", data = op)
vc_inc2_rob4 <- glm(y2 ~ okg + income + age + female + factor(educ) + factor(class5) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), family="binomial", data = op)

vc_inc2 <- broom::tidy(vc_inc2,conf.int=TRUE, exponentiate=FALSE)
vc_inc2 <- dplyr::filter(vc_inc2, term=="okg")
vc_inc2$model <- "main"

vc_inc2_rob1 <- broom::tidy(vc_inc2_rob1,conf.int=TRUE, exponentiate=FALSE)
vc_inc2_rob1 <- dplyr::filter(vc_inc2_rob1, term=="okg")
vc_inc2_rob1$model <- "main + fiscal constraint"

vc_inc2_rob2 <- broom::tidy(vc_inc2_rob2,conf.int=TRUE, exponentiate=FALSE)
vc_inc2_rob2 <- dplyr::filter(vc_inc2_rob2, term=="okg")
vc_inc2_rob2$model <- "main + unemp benefit"

vc_inc2_rob3 <- broom::tidy(vc_inc2_rob3,conf.int=TRUE, exponentiate=FALSE)
vc_inc2_rob3 <- dplyr::filter(vc_inc2_rob3, term=="okg")
vc_inc2_rob3$model <- "main + vulnerability"

vc_inc2_rob4 <- broom::tidy(vc_inc2_rob4,conf.int=TRUE, exponentiate=FALSE)
vc_inc2_rob4 <- dplyr::filter(vc_inc2_rob4, term=="okg")
vc_inc2_rob4$model <- "main + all 3"

vc_inc2_robust <- rbind(vc_inc2, vc_inc2_rob1, vc_inc2_rob2, vc_inc2_rob3, vc_inc2_rob4)

ref <- c("okg", 0, 0, 0, 0, 0, 0, "main")

vc_inc2_robust <- rbind(ref, vc_inc2_robust)

vc_inc2_robust$label <- NA
vc_inc2_robust$label[vc_inc2_robust$conf.high==0] <- "Incumbent (Reference)"
vc_inc2_robust$label[vc_inc2_robust$conf.high!=0] <- "Opposition"


vc_inc2_robust[,c(2:7)]  <-  apply(vc_inc2_robust[,c(2:7)], 2, function(x) as.numeric(as.character(x)));

ggplot(vc_inc2_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Kids and Incumbency") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()

# joint plot

fig1_controls <- ggplot(fig3a_mn_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Respondent") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted") +
  scale_color_grey()

fig2_controls <- ggplot(fig3b_mn_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Kids") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()

fig3_controls <- ggplot(vc_inc_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Respondent") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()

fig4_controls <- ggplot(vc_inc2_robust, aes(x=estimate,y=label, group=model, color=model, shape=model))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width=0.2)) + 
  xlab("Opportunity Kids") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")+
  scale_color_grey()

ggarrange(fig1_controls, fig2_controls, fig3_controls, fig4_controls, ncol=2, nrow=2, common.legend = TRUE, legend="bottom")
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D1.pdf"), width=9, height=6)

# Table D.1 ----

fig4_inc <- glm(y_inc~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_opp <- glm(y_opp~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_rad <- glm(y_rad~quadrant+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")


fig4_inc_robust <- glm(y_inc~quadrant+age_grp+female+educ_grp2+factor(class8) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), data=op, family="binomial")
fig4_opp_robust <- glm(y_opp~quadrant+age_grp+female+educ_grp2+factor(class8) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), data=op, family="binomial")
fig4_rad_robust <- glm(y_rad~quadrant+age_grp+female+educ_grp2+factor(class8) + factor(fiscal_constraint) + factor(ever_unemp_benefits) + factor(empstat_vuln) + factor(country), data=op, family="binomial")

# Table out


texreg(list(fig4_inc, fig4_inc_robust, fig4_opp, fig4_opp_robust, fig4_rad, fig4_rad_robust), digits=3, 
       label="tab:tab_controls_joint", caption="Opportunity Types and Party Support", caption.above=TRUE,
       custom.model.names=c("Inc", "Inc", "MS Opp",  "MS Opp", "Radical", "Radical"), 
       custom.coef.map = list('quadrantapprehensive' = 'Type: Apprehensive',
                              'quadrantaspirational' = 'Type: Aspirational',
                              'quadrantburdened' = 'Type: Burdened',
                              'age_grp26-35' = 'Age: 26-35',
                              'age_grp36-45' = 'Age: 36-45',
                              'age_grp46-55' = 'Age: 46-55',
                              'age_grp56-65' = 'Age: 56-65',
                              'age_grp66+'='Age: 66 and older',
                              'female1' = 'Female',
                              'educ_grp2.L' = 'Educ: Medium',
                              'educ_grp2.Q' = 'Educ: High',
                              'factor(class8)employers' = 'Class: Employers',
                              'factor(class8)managers' = 'Class: Managers',
                              'factor(class8)prod workers' = 'Class: Production Workers',
                              'factor(class8)service workers' = 'Class: Service Workers',
                              'factor(class8)small bus owners' = 'Class: Small Business Owners',
                              'factor(class8)sociocult profs' = 'Class: Soc. Cult. Profs',
                              'factor(class8)technical profs' = 'Class: Tech. Profs',
                              'factor(fiscal_constraint)1' = 'Fiscal Constraint: Agree',
                              'factor(ever_unemp_benefits)1' = 'Ever received unemp benefits', 
                              'factor(empstat_vuln)1' = 'Current LM Vulnerability'
                              ),
       custom.note=paste("%stars. All models include country-fixed effects."),
       file=paste0(outpath, "BJPS_Appendix_Table_D1.tex"))

# Figure D.2 ----

op$yg <- as.factor(op$yg)
op$yg <- relevel(op$yg, ref="1 inc")

# Opportunity Respondent

fig3a_green_mn1_mod <- multinom(yg ~ org + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

fig3a_green_mn1 <- broom::tidy(fig3a_green_mn1_mod,conf.int=TRUE, exponentiate=FALSE) # countries without green parties create warnings
fig3a_green_mn1 <- dplyr::filter(fig3a_green_mn1, term=="org")

ref3a_green_mn1 <- c("1 inc", "org", 0, 0, 0, 0, 0, 0)
fig3a_green_label <- c("Mainstream Opposition", "Green", "Radical", "Incumbent (Reference)")

fig3a_green_plot1 <- rbind(fig3a_green_mn1, ref3a_green_mn1)
fig3a_green_plot1 <- cbind(fig3a_green_label, fig3a_green_plot1)

fig3a_green_plot1 %<>% arrange(fig3a_green_label)
fig3a_green_plot1[,c(4:9)]  <-  apply(fig3a_green_plot1[,c(4:9)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3a_green_plot1, aes(x=estimate,y=fig3a_green_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
   xlim((-0.25), 0.25) +
  xlab("Opportunity Respondent and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D2a.pdf"), width=6, height=4)

# Opportunity Kids

fig3a_green_mn2_mod <- multinom(yg ~ okg + income + age + female + factor(educ) + factor(class5) + factor(country), data = op)

fig3a_green_mn2 <- broom::tidy(fig3a_green_mn2_mod,conf.int=TRUE, exponentiate=FALSE) # some issues with countries without green parties..
fig3a_green_mn2 <- dplyr::filter(fig3a_green_mn2, term=="okg")

ref3a_green_mn2 <- c("1 inc", "okg", 0, 0, 0, 0, 0, 0)

fig3a_green_plot2 <- rbind(fig3a_green_mn2, ref3a_green_mn2)
fig3a_green_plot2 <- cbind(fig3a_green_label, fig3a_green_plot2)

fig3a_green_plot2 %<>% arrange(fig3a_green_label)
fig3a_green_plot2[,c(4:9)]  <-  apply(fig3a_green_plot2[,c(4:9)], 2, function(x) as.numeric(as.character(x)));

ggplot(fig3a_green_plot2, aes(x=estimate,y=fig3a_green_label))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + 
   xlim((-0.25), 0.25) +
  xlab("Opportunity Kids and Party Support") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D2b.pdf"), width=6, height=4)

# Figure D.3 ----

# define emblematic cases

typ_comf <- output %>% dplyr::filter(quadrant=="comfortable", age_grp=="56-65", educ_grp2=="high", female=="0", class8=="managers") %>%
  group_by(party) %>% summarise(mean=mean(mean))
typ_comf$quadrant <- "56-65yo high-skill male manager"

typ_app <- output %>% dplyr::filter(quadrant=="apprehensive", age_grp=="46-55", educ_grp2=="medium", female=="0", class8=="clerks") %>%
  group_by(party) %>% summarise(mean=mean(mean))
typ_app$quadrant <- "46-55yo mid-skill male clerk"

typ_asp <- output %>% dplyr::filter(quadrant=="aspirational", age_grp=="18-25", educ_grp2=="high", female=="1", class8=="service workers") %>%
  group_by(party) %>% summarise(mean=mean(mean))
typ_asp$quadrant <- "18-25yo high-skill female service"

typ_bur <- output %>% dplyr::filter(quadrant=="burdened", age_grp=="26-35", educ_grp2=="low", female=="0", class8=="service workers") %>%
  group_by(party) %>% summarise(mean=mean(mean))
typ_bur$quadrant <- "26-35yo low-skill male service"

typical <- rbind(typ_comf, typ_app, typ_asp, typ_bur)

bar <- output_avg %>% dplyr::select(party, bar)
bar <- unique(bar)

typical <- merge(typical, bar, by="party")

typical$quadrant <-  ordered(typical$quadrant, levels = c("56-65yo high-skill male manager", "46-55yo mid-skill male clerk", "18-25yo high-skill female service", "26-35yo low-skill male service"))

ggplot(typical, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=5) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant) + ylim(0,0.55) +
  scale_shape_manual(name="", labels = c("Incumbent", "Mainstream Opposition", "Radical"), values=c(10,4,13)) +
  labs(x="",
       y="",
       title="Predicted probabilities of emblematic cases") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Appendix_Figure_D3.pdf"), height=6, width=8)

# Figure D.4 ----

op$yy <- as.factor(op$y)
op$yy <- relevel(op$yy, ref="2 opp")

op$qy <- relevel(op$quadrant, ref="aspirational")

fig4_mn <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op)

fig4_mn_plot <- broom::tidy(fig4_mn,conf.int=TRUE, exponentiate=FALSE)
fig4_mn_plot <- dplyr::filter(fig4_mn_plot, term%in%c("qycomfortable", "qyapprehensive", "qyburdened"))

fig4_mn_plot$facet <- ifelse(fig4_mn_plot$y.level=="1 inc", "Incumbent", "Radical")
fig4_mn_plot$quadrant <- NA
fig4_mn_plot$quadrant[fig4_mn_plot$term=="qycomfortable"] <- "Comfortable"
fig4_mn_plot$quadrant[fig4_mn_plot$term=="qyapprehensive"] <- "Apprehensive"
fig4_mn_plot$quadrant[fig4_mn_plot$term=="qyburdened"] <- "Burdened"

ggplot(fig4_mn_plot, aes(x=estimate,y=quadrant))+
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high)) + facet_wrap(~facet) +
  xlab("Relative to (1) aspirational (2) opposition voters") + ylab("") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D4.pdf"), height=4, width=6)

# Figure D.5 ----

fig4_mn_DK <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Denmark"))
fig4_mn_DE <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Germany"))
fig4_mn_IE <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Ireland"))
fig4_mn_IT <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Italy"))
fig4_mn_NL <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Netherlands"))
fig4_mn_ES <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Spain"))
fig4_mn_SE <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="Sweden"))
fig4_mn_UK <- multinom(yy~qy+age_grp+female+educ_grp2+factor(class8), data=subset(op, country!="United Kingdom"))

fig4_mn_DK <- broom::tidy(fig4_mn_DK,conf.int=TRUE)
fig4_mn_DK$cntry <- "excl. DK"
fig4_mn_DE <- broom::tidy(fig4_mn_DE,conf.int=TRUE)
fig4_mn_DE$cntry <- "excl. DE"
fig4_mn_IE <- broom::tidy(fig4_mn_IE,conf.int=TRUE)
fig4_mn_IE$cntry <- "excl. IE"
fig4_mn_IT <- broom::tidy(fig4_mn_IT,conf.int=TRUE)
fig4_mn_IT$cntry <- "excl. IT"
fig4_mn_NL <- broom::tidy(fig4_mn_NL,conf.int=TRUE)
fig4_mn_NL$cntry <- "excl. NL"
fig4_mn_ES <- broom::tidy(fig4_mn_ES,conf.int=TRUE)
fig4_mn_ES$cntry <- "excl. ES"
fig4_mn_SE <- broom::tidy(fig4_mn_SE,conf.int=TRUE)
fig4_mn_SE$cntry <- "excl. SE"
fig4_mn_UK <- broom::tidy(fig4_mn_UK,conf.int=TRUE)
fig4_mn_UK$cntry <- "excl. UK"

fig4_mn_cntry <- do.call("rbind", list(fig4_mn_DK, fig4_mn_DE, fig4_mn_IE, fig4_mn_IT, fig4_mn_NL, fig4_mn_ES, fig4_mn_SE, fig4_mn_UK))

fig4_mn_cntry_plot <- dplyr::filter(fig4_mn_cntry, term%in%c("qycomfortable", "qyapprehensive", "qyburdened"))

fig4_mn_cntry_plot$facet <- ifelse(fig4_mn_cntry_plot$y.level=="1 inc", "Incumbent", "Radical")
fig4_mn_cntry_plot$quadrant <- NA
fig4_mn_cntry_plot$quadrant[fig4_mn_cntry_plot$term=="qycomfortable"] <- "Comfortable"
fig4_mn_cntry_plot$quadrant[fig4_mn_cntry_plot$term=="qyapprehensive"] <- "Apprehensive"
fig4_mn_cntry_plot$quadrant[fig4_mn_cntry_plot$term=="qyburdened"] <- "Burdened"

ggplot(fig4_mn_cntry_plot, aes(x=estimate,y=quadrant, color=factor(cntry), group=factor(cntry))) +
  geom_pointrange(aes(xmin=conf.low,
                     xmax=conf.high), position=position_dodge(width = .3)) + facet_wrap(~facet) +
  xlab("Relative to (1) aspirational (2) opposition voters") + ylab("") + scale_color_discrete(name = "") +
  geom_vline(aes(xintercept=0), color="red", linetype="dotted")
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D5.pdf"), height=4, width=6)

# Figure D.6 ----

op$female <- as.factor(op$female)
op$age_grp <- as.factor(op$age_grp)
op$quadrant3 <- as.factor(op$quadrant3)

fig4_q3_inc <- glm(y_inc~quadrant3+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_q3_opp <- glm(y_opp~quadrant3+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_q3_rad <- glm(y_rad~quadrant3+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Calculate average predicted probabilities

output_q3_inc = glm.predict::predicts(fig4_q3_inc, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_q3_inc$party <- "Incumbent"
output_q3_opp = predicts(fig4_q3_opp, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_q3_opp$party <- "Mainstream Opposition"
output_q3_rad = predicts(fig4_q3_rad, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_q3_rad$party <- "Radical"

output_q3<- full_join(output_q3_inc, output_q3_opp)
output_q3<- full_join(output_q3, output_q3_rad)

output_q3<-  output_q3 %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_q3_avg <- output_q3 %>% group_by(party, quadrant3) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output_q3_avg$quadrant3 <-  ordered(output_q3_avg$quadrant3, levels = c("highincpos", "highincneg", "medincpos", "medincneg", "lowincpos", "lowincneg"))

ggplot(output_q3_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=5) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant3, nrow=3)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Incumbent", "Mainstream Opposition", "Radical")), values=c(10,4,13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())
ggsave(paste0(outpath, "BJPS_Appendix_Figure_D6.pdf"), height=6, width=8)

# Figure D.7 ----

fig4_kids_inc <- glm(y_inc~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_kids_opp <- glm(y_opp~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_kids_rad <- glm(y_rad~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Calculate average predicted probabilities

output_kids_inc = predicts(fig4_kids_inc, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_inc$party <- "incumbent"
output_kids_opp = predicts(fig4_kids_opp, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_opp$party <- "Mainstream Opposition"
output_kids_rad = predicts(fig4_kids_rad, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_rad$party <- "Radical"

output_kids <- full_join(output_kids_inc, output_kids_opp)
output_kids <- full_join(output_kids, output_kids_rad)

output_kids <-  output_kids %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_kids_avg <- output_kids %>% group_by(party, quadrant_kids) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output_kids_avg$quadrant_kids <-  ordered(output_kids_avg$quadrant_kids, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

ggplot(output_kids_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=5) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant_kids)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Incumbent", "Mainstream Opposition", "Radical")), values=c(10,4,13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Appendix_Figure_D7.pdf"), height=6, width=8)

# Figure D.8 ----

op$female <- as.factor(op$female)
op$age_grp <- as.factor(op$age_grp)

fig4_lm_inc <- glm(y_inc~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_lm_opp <- glm(y_opp~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_lm_rad <- glm(y_rad~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Calculate average predicted probabilities

output_lm_inc = predicts(fig4_lm_inc, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_inc$party <- "incumbent"
output_lm_opp = predicts(fig4_lm_opp, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_opp$party <- "Mainstream Opposition"
output_lm_rad = predicts(fig4_lm_rad, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_rad$party <- "Radical"

output_lm <- full_join(output_lm_inc, output_lm_opp)
output_lm <- full_join(output_lm, output_lm_rad)

output_lm <-  output_lm %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_lm_avg <- output_lm %>% group_by(party, quadrant_lm) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output_lm_avg$quadrant_lm <-  ordered(output_lm_avg$quadrant_lm, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))

ggplot(output_lm_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=5) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant_lm)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Incumbent", "Mainstream Opposition", "Radical")), values=c(10,4,13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Appendix_Figure_D8.pdf"), height=6, width=8)

# Figure D.9 ----

fig4_kids_rr <- glm(rr~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_kids_msr <- glm(msr~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_kids_msl <- glm(msl~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_kids_rl <- glm(rl~quadrant_kids+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Calculate average predicted probabilities

output_kids_rr = predicts(fig4_kids_rr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_rr$party <- "Rad Right"
output_kids_msr = predicts(fig4_kids_msr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_msr$party <- "MS Right"
output_kids_msl = predicts(fig4_kids_msl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_msl$party <- "MS Left"
output_kids_rl = predicts(fig4_kids_rl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_kids_rl$party <- "Rad Left"

output_kids2 <- full_join(output_kids_rr, output_kids_msr)
output_kids2 <- full_join(output_kids2, output_kids_msl)
output_kids2 <- full_join(output_kids2, output_kids_rl)


output_kids2 <-  output_kids2 %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_kids2_avg <- output_kids2 %>% group_by(party, quadrant_kids) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output_kids2_avg$quadrant_kids <-  ordered(output_kids2_avg$quadrant_kids, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))
output_kids2_avg$party <-  ordered(output_kids2_avg$party, levels = c("Rad Left", "MS Left", "MS Right", "Rad Right"))

ggplot(output_kids2_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=6) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant_kids)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Rad Left", "MS Left", "MS Right", "Rad Right")), values=c(10,4, 8, 13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Appendix_Figure_D9.pdf"), width=8, height=6)

# Figure D.10 ----

fig4_lm_rr <- glm(rr~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_lm_msr <- glm(msr~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_lm_msl <- glm(msl~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")
fig4_lm_rl <- glm(rl~quadrant_lm+age_grp+female+educ_grp2+factor(class8)+factor(country), data=op, family="binomial")

# Calculate average predicted probabilities

output_lm_rr = predicts(fig4_lm_rr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_rr$party <- "Rad Right"
output_lm_msr = predicts(fig4_lm_msr, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_msr$party <- "MS Right"
output_lm_msl = predicts(fig4_lm_msl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_msl$party <- "MS Left"
output_lm_rl = predicts(fig4_lm_rl, "F;F;F;F;F;F", sim.count = 10000, set.seed = 123)
output_lm_rl$party <- "Rad Left"

output_lm2 <- full_join(output_lm_rr, output_lm_msr)
output_lm2 <- full_join(output_lm2, output_lm_msl)
output_lm2 <- full_join(output_lm2, output_lm_rl)


output_lm2 <-  output_lm2 %>% rename(class8=`factor(class8`,
                             country=`factor(country`)

output_lm2_avg <- output_lm2 %>% group_by(party, quadrant_lm) %>% # unweighted average predicted probabilities
  summarise(mean=mean(mean)) %>% ungroup() %>%
  group_by(party) %>%
  mutate(bar=mean(mean))

# Plot average predicted probabilities

output_lm2_avg$quadrant_lm <-  ordered(output_lm2_avg$quadrant_lm, levels = c("comfortable", "apprehensive", "aspirational", "burdened"))
output_lm2_avg$party <-  ordered(output_lm2_avg$party, levels = c("Rad Left", "MS Left", "MS Right", "Rad Right"))

ggplot(output_lm2_avg, aes(x=party)) + 
  geom_bar(aes(y=bar), color="black", fill="white", stat="identity") + 
  geom_point(aes(y=mean, shape=party), size=6) + 
  geom_segment(aes(x=party, xend=party, y=bar, yend=mean)) +
  facet_wrap(~quadrant_lm)  + ylim(0,0.55) +
  scale_shape_manual(name="", labels = (c("Rad Left", "MS Left", "MS Right", "Rad Right")), values=c(10,4, 8, 13)) +
  labs(x="",
       y="",
       title="Average predicted probabilities") +
  theme(legend.position = "bottom", axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggsave(paste0(outpath, "BJPS_Appendix_Figure_D10.pdf"), width=8, height=6)

# Figure E.1 ----

w_ms <- op %>% filter(!is.na(quadrant)) %>%
  filter(y_rad==0) %>%
                     mutate(one=1, wone=Weight_Age_sex___education) %>%
  group_by(country) %>% mutate(cntrytotal=sum(one), wcntrytotal=sum(Weight_Age_sex___education, na.rm=T)) %>%
  group_by(country, quadrant) %>%
  summarise(total=sum(one, na.rm=T), wtotal=sum(wone, na.rm=T), 
         share=total/cntrytotal, wshare=wtotal/wcntrytotal,
         perc=round(share*100, 1), wperc=round(wshare*100, 1)) %>%
  filter(row_number(total) == 1) %>%
    ungroup()

w_ms_table <- w_ms %>% dplyr::select(country, quadrant, wperc)

w_ms_table %>% filter(quadrant=="apprehensive" | quadrant=="burdened") %>%
  group_by(country) %>% summarise(threat=sum(wperc)) %>%
  ggplot(aes(x=reorder(country, threat), y=threat)) + geom_bar(stat="identity")

w_ms_table %>% filter(quadrant=="apprehensive" | quadrant=="burdened") %>%
  group_by(country) %>%
  ggplot(aes(x=reorder(country, wperc), y=wperc, group=quadrant, fill=quadrant)) + geom_bar(stat="identity", position="stack") +
  xlab("") + ylab("Relative Share among Mainstream Party Voters") + scale_fill_grey(name="")

ggsave(paste0(outpath, "BJPS_Appendix_Figure_E1.pdf"), height=6, width=9)

# ENDEND.
